
rm(list=ls())



pkgs <- c("dplyr", "ggplot2", "cjoint", "cregg")

# Load all above packages. Install if they have not been installed.
usePackage <- function(p){
  for (pkg in p){
    if (!is.element(pkg, installed.packages()[,1]))
      install.packages(pkg, dep = TRUE, repos = "https://cloud.r-project.org/")
    require(pkg, character.only = TRUE)
  }
}
usePackage(pkgs)

## Must set working directory  

## Load the data
d0 <- read.csv("CPConjointData.csv")

## Code the conjoint tasks 

table(d0$Conjoint1_SC, useNA = "always")
table(d0$Conjoint2_SC, useNA = "always")
table(d0$Conjoint3_SC, useNA = "always")
d0$Conjoint1_SC <- ifelse(d0$Conjoint1_SC == "Scenario A", 1, 0)
d0$Conjoint2_SC <- ifelse(d0$Conjoint2_SC == "Scenario A", 1, 0)
d0$Conjoint3_SC <- ifelse(d0$Conjoint3_SC == "Scenario A", 1, 0)



## Code partisanship and education

d0$Republican <- 0
d0$Republican[d0$Demo_party_SC == "Republican " | d0$Demo_partylean_SC == "Republican Party"] <- 1
d0$Democrat <- 0
d0$Democrat[d0$Demo_party_SC == "Democrat " | d0$Demo_partylean_SC == "Democratic Party "] <- 1

d0$Independent <- as.numeric(ifelse(d0$Democrat == 1 | d0$Republican == 1, 0, 1))
sum(d0$Independent)

d0$CollegeDegree <- ifelse(d0$Demo_education_SC == "College graduate" | d0$Demo_education_SC == "Graduate degree", "Yes", "No")


d1 <- subset(d0, d0$Independent !=1)

d1$InParty <- d1$Inparty2_QD


### Code conjoint trials:

### National election

# Trial 1, Scenario A
d1$Nat_elect_QD_1A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_1A)
d1$Nat_elect_QD_1A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_1A)

# Trial 1, Scenario B
d1$Nat_elect_QD_1B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_1B)
d1$Nat_elect_QD_1B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_1B)

# Trial 2, Scenario A
d1$Nat_elect_QD_2A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_2A)
d1$Nat_elect_QD_2A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_2A)

# Trial 2, Scenario B
d1$Nat_elect_QD_2B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_2B)
d1$Nat_elect_QD_2B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_2B)

# Trial 3, Scenario A
d1$Nat_elect_QD_3A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_3A)
d1$Nat_elect_QD_3A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_3A)

# Trial 3, Scenario B
d1$Nat_elect_QD_3B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Nat_elect_QD_3B)
d1$Nat_elect_QD_3B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Nat_elect_QD_3B)

### Local election

# Trial 1, Scenario A
d1$Local_elect_QD_1A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_1A)
d1$Local_elect_QD_1A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_1A)

# Trial 1, Scenario B
d1$Local_elect_QD_1B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_1B)
d1$Local_elect_QD_1B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_1B)

# Trial 2, Scenario A
d1$Local_elect_QD_2A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_2A)
d1$Local_elect_QD_2A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_2A)

# Trial 2, Scenario B
d1$Local_elect_QD_2B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_2B)
d1$Local_elect_QD_2B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_2B)

# Trial 3, Scenario A
d1$Local_elect_QD_3A <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_3A)
d1$Local_elect_QD_3A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_3A)

# Trial 3, Scenario B
d1$Local_elect_QD_3B <- gsub("\\$\\{e://Field/Inparty_QD\\}", "In-party", d1$Local_elect_QD_3B)
d1$Local_elect_QD_3B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Local_elect_QD_3B)

### Election Legislation

# Trial 1, Scenario A
d1$Leg_change_QD_1A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_1A)
d1$Leg_change_QD_1A <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_1A)

# Trial 1, Scenario B
d1$Leg_change_QD_1B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_1B)
d1$Leg_change_QD_1B <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_1B)

# Trial 2, Scenario A
d1$Leg_change_QD_2A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_2A)
d1$Leg_change_QD_2A <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_2A)

# Trial 2, Scenario B
d1$Leg_change_QD_2B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_2B)
d1$Leg_change_QD_2B <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_2B)

# Trial 3, Scenario A
d1$Leg_change_QD_3A <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_3A)
d1$Leg_change_QD_3A <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_3A)

# Trial 3, Scenario B
d1$Leg_change_QD_3B <- gsub("\\$\\{e://Field/Outparty_QD\\}", "Out-party", d1$Leg_change_QD_3B)
d1$Leg_change_QD_3B <- gsub("\\$\\{e://Field/Outparty2_QD\\}", "Out-party", d1$Leg_change_QD_3B)

d <- dplyr::mutate(d1, ID = dplyr::row_number())

################
list1 <- list(
  NatElect = list(
    ~ Nat_elect_QD_1A + Nat_elect_QD_2A + Nat_elect_QD_3A ,
    ~ Nat_elect_QD_1B + Nat_elect_QD_2B + Nat_elect_QD_3B ),
  LocalElect = list(
    ~ Local_elect_QD_1A + Local_elect_QD_2A + Local_elect_QD_3A ,
    ~ Local_elect_QD_1B + Local_elect_QD_2B + Local_elect_QD_3B),
  LegChange = list(
    ~ Leg_change_QD_1A+ Leg_change_QD_2A +Leg_change_QD_3A ,
    ~ Leg_change_QD_1B+ Leg_change_QD_2B+Leg_change_QD_3B),
  SalChange = list(
    ~ Sal_change_QD_1A+Sal_change_QD_2A+Sal_change_QD_3A,
    ~ Sal_change_QD_1B+Sal_change_QD_2B+Sal_change_QD_3B)
)
# task variables
list2 <- list(choice = ~ Conjoint1_SC + Conjoint2_SC + Conjoint3_SC )


str(long <- cregg::cj_tidy(d,
                    profile_variables = list1,
                    task_variables = list2,
                    id = ~ ID))
stopifnot(nrow(long) == nrow(d)*3*2)
dim(long)
table(long$profile, useNA = "always")
table(long$choice, useNA = "always")

long$chosen <- ifelse((long$profile == "A" & long$choice == 1) | 
                        (long$profile == "B" & long$choice == 0), 1, 0)
table(long$chosen, useNA = "always")
dim(long)
long <- subset(long, long$chosen != "NA")
dim(long)

long$NatElect <- gsub("Out-party win", "Out-party wins", long$NatElect)
long$NatElect <- gsub("In-party win", "In-party wins", long$NatElect)
long$NatElect <- gsub("take control of Congress", "takes control of Congress", long$NatElect)
long$LocalElect <- gsub("take control of", "takes control of", long$LocalElect)
long$LegChange <- gsub("a Out-party candidate", "an out-party candidate", long$LegChange)
long$LegChange <- gsub("all Out-party", "all out-party members", long$LegChange)
long$LegChange <- gsub("some Out-party", "some out-party members", long$LegChange)
long$LegChange <- gsub("restricts a Out-party candidate", "restricts an out-party candidate", long$LegChange)

table(long$NatElect, useNA = "always")
long$Presidential_Election_Outcome_2024 <- factor(long$NatElect, levels = c("Out-party wins the Presidency and takes control of Congress.",
                                                                            "Out-party wins the Presidency but Congress is divided.",
                                                                            "In-party wins the Presidency but Congress is divided.",
                                                                            "In-party wins the Presidency and takes control of Congress."))
table(long$Presidential_Election_Outcome_2024, useNA = "always")

long$Local_Election_Outcome_2024 <- factor(long$LocalElect, levels = c("Out-party takes control of your local governing board.",
                                                                       "Your local governing board is evenly split between Democrats and Republicans.",
                                                                       "In-party takes control of your local governing board."))
table(long$Local_Election_Outcome_2024, useNA = "always")

long$Election_Legislation <- factor(long$LegChange, levels = c("The state legislature recommends a cap on the number of campaign volunteers, applying to all candidates in state and local campaigns.",
                                                               "The state legislature restricts an out-party candidate from advertising in state and local campaigns in response to recent controversies.",
                                                               "The state legislature bans an out-party candidate from competing in state and local campaigns in response to recent controversies.",
                                                               "The state legislature bans some out-party members from competing in state and local campaigns in response to recent controversies.",
                                                               "The state legislature bans all out-party members from competing in state and local campaigns in response to recent controversies."))
table(long$Election_Legislation, useNA = "always")

long$Salary_Change <- factor(long$SalChange, levels = c("The state legislature leaves your pay unchanged.",
                                                        "The state legislature cuts your pay by 20%.",
                                                        "The state legislature cuts your pay by 10%.",
                                                        "The state legislature cuts your pay by 5%.",
                                                        "The state legislature raises your pay by 5%.",
                                                        "The state legislature raises your pay by 10%.",
                                                        "The state legislature raises your pay by 20%."))
table(long$Salary_Change, useNA = "always")

long$ID <- as.factor(long$ID)
table(long$ID)
table(table(long$ID))

#####order factors######

long <- within(long, Presidential_Election_Outcome_2024 <- relevel(Presidential_Election_Outcome_2024, ref = "Out-party wins the Presidency and takes control of Congress."))  
long <- within(long, Local_Election_Outcome_2024 <- relevel(Local_Election_Outcome_2024, ref = "Your local governing board is evenly split between Democrats and Republicans."))
long <- within(long, Election_Legislation <- relevel(Election_Legislation, ref = "The state legislature recommends a cap on the number of campaign volunteers, applying to all candidates in state and local campaigns."))
long <- within(long, Salary_Change <- relevel(Salary_Change, ref = "The state legislature leaves your pay unchanged."))

## AMCE
detach("package:cregg", unload = TRUE)

f1Main <- cjoint::amce(chosen ~ Presidential_Election_Outcome_2024 + Local_Election_Outcome_2024 + Election_Legislation + Salary_Change, data = long, cluster = TRUE, respondent.id = "ID")
summary(f1Main)

## Baseline Selection Probability

long$chosen <- as.numeric(long$chosen)


baseline_subset <- long %>%
  filter(Election_Legislation == "The state legislature bans all out-party members from competing in state and local campaigns in response to recent controversies." &
           Presidential_Election_Outcome_2024 == "In-party wins the Presidency and takes control of Congress.")

baseline_prob <- mean(baseline_subset$chosen, na.rm = TRUE)

print(baseline_prob)

## Using weights 

f1MainWeighted <- cjoint::amce(chosen ~Presidential_Election_Outcome_2024 + Local_Election_Outcome_2024 + Election_Legislation  + Salary_Change, data= long, cluster=TRUE, respondent.id= "ID", weights = "Weight")
summary(f1MainWeighted)

## Plot Figure 1

attributes <- c(
  "Election Legislation", "Election Legislation", "Election Legislation", "Election Legislation",
  "2024 Presidential Election Outcome", "2024 Presidential Election Outcome", "2024 Presidential Election Outcome",
  "2024 Local Election Outcome", "2024 Local Election Outcome",
  "Salary Change", "Salary Change", "Salary Change", "Salary Change", "Salary Change", "Salary Change"
)

levels <- c(
  "Restricts out-party candidate advertising", "Bans out-party candidate", "Bans some out-party members", "Bans all out-party members",
  "Out-party wins Presidency, divided Congress", "In-party wins Presidency, divided Congress", "In-party wins Presidency & Congress", # This list of categories is wrong unless we change the reference category earlier in the script. I have now made that change earlier in the script.
  "Out-party takes local board", "In-party takes local board",
  "Legislature cuts pay by 20%", "Legislature cuts pay by 10%", "Legislature cuts pay by 5%",
  "Legislature raises pay by 5%", "Legislature raises pay by 10%", "Legislature raises pay by 20%"
)

summary(f1Main)

# These come from summary(f1Main)
estimates <- c(-0.0437954, -0.0787197, -0.1217270, -0.2022350,
               0.2244558, 0.4765559, 0.5673831,  
               -0.1558702, -0.0186192,
               -0.0627786, -0.0195846, -0.0070590, 0.0188375, -0.0051269, -0.0468885)

std_errs <- c(0.031198, 0.031386, 0.031219, 0.032309,
              0.025612, 0.028414, 0.027355, 
              0.025144, 0.024566,
              0.034115, 0.034443, 0.034459, 0.032631, 0.035166, 0.031558)


# Calculate the 95% confidence intervals
lower_bounds <- estimates - (1.96 * std_errs)
upper_bounds <- estimates + (1.96 * std_errs)

# Combine into a data frame
results_df <- data.frame(
  Attribute = attributes,
  Level = levels,
  Estimate = estimates,
  StdErr = std_errs,
  LowerCI = lower_bounds,
  UpperCI = upper_bounds
)

# Add 'Significant' column based on the 95% confidence interval
results_df$Significant <- with(results_df, LowerCI > 0 | UpperCI < 0)

# Update baseline levels labels
baseline_levels <- data.frame(
  Attribute = c("Election Legislation", "2024 Presidential Election Outcome", 
                "2024 Local Election Outcome", "Salary Change"),
  Level = c("Non-binding cap on campaign volunteers", "Out-party wins Presidency & Congress", 
            "Local board evenly divided", "Legislature leaves pay unchanged"),
  Estimate =0,
  StdErr = 0,
  LowerCI = 0,
  UpperCI = 0,
  Significant = FALSE   
)


baseline_levels$Level <- gsub("Non-binding cap on campaign volunteers", "Non-binding cap on\ncampaign volunteers", baseline_levels$Level)
baseline_levels$Level <- gsub("Out-party wins Presidency & Congress", "Out-party wins\nPresidency & Congress", baseline_levels$Level)
baseline_levels$Level <- gsub("Local board evenly divided", "Local board\nevenly divided", baseline_levels$Level)
baseline_levels$Level <- gsub("Legislature leaves pay unchanged", "Pay unchanged", baseline_levels$Level)

results_with_baseline <- rbind(results_df, baseline_levels)


# The order of the attributes given here determines the order in which they appear in the plot:
results_with_baseline$Attribute <- factor(results_with_baseline$Attribute, 
                                          levels = c("Election Legislation", "2024 Presidential Election Outcome", "2024 Local Election Outcome", "Salary Change"))
table(results_with_baseline$Attribute, useNA = "always")


# Add line breaks to the levels for better visualization
results_with_baseline$Level <- gsub("Restricts out-party candidate advertising", "Restricts out-party\ncandidate advertising", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Bans out-party candidate", "Bans out-party\ncandidate", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Bans some out-party members", "Bans some\nout-party members", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Bans all out-party members", "Bans all\nout-party members", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Out-party wins Presidency, divided Congress", "Out-party wins Presidency,\ndivided Congress", results_with_baseline$Level)
results_with_baseline$Level <- gsub("In-party wins Presidency, divided Congress", "In-party wins Presidency,\ndivided Congress", results_with_baseline$Level)
results_with_baseline$Level <- gsub("In-party wins Presidency & Congress", "In-party wins Presidency\n& Congress", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Out-party takes local board", "Out-party takes\nlocal board", results_with_baseline$Level)
results_with_baseline$Level <- gsub("In-party takes local board", "In-party takes\nlocal board", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature cuts pay by 20%", "Pay cut by 20%", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature cuts pay by 10%", "Pay cut by 10%", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature cuts pay by 5%", "Pay cut by 5%", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature raises pay by 5%", "Pay raised by 5%", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature raises pay by 10%", "Pay raised by 10%", results_with_baseline$Level)
results_with_baseline$Level <- gsub("Legislature raises pay by 20%", "Pay raised by 20%", results_with_baseline$Level)


desired_order <- rev(c(
  "Non-binding cap on\ncampaign volunteers",
  "Restricts out-party\ncandidate advertising",
  "Bans out-party\ncandidate",
  "Bans some\nout-party members",
  "Bans all\nout-party members",
  "In-party wins Presidency\n& Congress",
  "In-party wins Presidency,\ndivided Congress",
  "Out-party wins Presidency,\ndivided Congress",
  "Out-party wins\nPresidency & Congress",
  "Local board\nevenly divided",
  "In-party takes\nlocal board",
  "Out-party takes\nlocal board",
  "Pay unchanged",
  "Pay cut by 20%",
  "Pay cut by 10%",
  "Pay cut by 5%",
  "Pay raised by 5%",
  "Pay raised by 10%",
  "Pay raised by 20%"
))
results_with_baseline$Level <- factor(results_with_baseline$Level, levels = desired_order)


ggplot(results_with_baseline, aes(x = Level, y = Estimate, fill = Significant)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("grey", NA)) +
  coord_flip() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.4, color = "black") +
  facet_wrap(~Attribute, scales = "free_y", ncol = 1) +
  labs(x = "", y = "AMCE Estimate", title = "") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 12, color = "black"),          
    axis.text.y = element_text(size = 12, color = "black"),          
    axis.title.x = element_blank(),
    axis.title.y = element_text(size = 14, color = "black"),         
    strip.text = element_text(size = 13, angle = 0, color = "black"), 
    panel.grid.major = element_blank(),             
    panel.grid.minor = element_blank(),             
    panel.background = element_rect(fill = "white", colour = "white"),
    axis.line.y = element_line(color = "black"),
    axis.ticks.y = element_line(color = "black"),
    axis.line = element_line(colour = "black"),
    legend.text = element_text(size = 12, color = "black"),          
    legend.title = element_text(size = 12, color = "black"),
    legend.position = "none"
  ) 
ggsave("Figure1.pdf", width = 6, height = 11, units = "in", dpi = 500)

## Marginal Means:

library(cregg)

# Rename columns in the data frame

names(long)[names(long) == "Local_Election_Outcome_2024"] <- "2024_Local_Election_Outcome"
names(long)[names(long) == "Presidential_Election_Outcome_2024"] <- "2024_Presidential_Election_Outcome"


mm_results <- mm(long, chosen ~ Election_Legislation + `2024_Local_Election_Outcome` + 
                   `2024_Presidential_Election_Outcome` + Salary_Change, id = ~ ID)

mm_data_frame <- mm_results

# Filter to keep only relevant columns
mm_data_frame <- mm_data_frame[, c("feature", "level", "estimate", "std.error", "lower", "upper")]

# Rename the columns for clarity
colnames(mm_data_frame) <- c("Attribute", "Level", "Estimate", "StdError", "LowerCI", "UpperCI")
mm_data_frame$Attribute <- gsub("_", " ", mm_data_frame$Attribute) # Replace underscores with spaces

long$Party <- as.factor(long$InParty)

mm_resultsPID <- cj(long, chosen ~ Election_Legislation + `2024_Local_Election_Outcome` + `2024_Presidential_Election_Outcome` + Salary_Change, estimate = "mm", id = ~ ID, by = ~ Party)

mm_resultsPID_df <- as.data.frame(mm_resultsPID)

mm_resultsPID_df <- mm_resultsPID_df[, c("feature", "level", "estimate", "std.error", "lower", "upper", "Party")]

# Rename the columns
colnames(mm_resultsPID_df) <- c("Attribute", "Level", "Estimate", "StdError", "LowerCI", "UpperCI", "Party")
mm_resultsPID_df$Attribute <- gsub("_", " ", mm_resultsPID_df$Attribute) 

mm_resultsPID_df$Attribute <- factor(mm_resultsPID_df$Attribute, 
                                     levels = c("Election Legislation", "2024 Presidential Election Outcome", "2024 Local Election Outcome", "Salary Change"))


new_labels <- rev(c(
  "The state legislature recommends a cap on the number of campaign volunteers, applying to all candidates in state and local campaigns." = "Non-binding cap on\ncampaign volunteers",
  "The state legislature restricts an out-party candidate from advertising in state and local campaigns in response to recent controversies." = "Restricts out-party\ncandidate advertising",
  "The state legislature bans an out-party candidate from competing in state and local campaigns in response to recent controversies." = "Bans out-party\ncandidate",
  "The state legislature bans some out-party members from competing in state and local campaigns in response to recent controversies." = "Bans some\nout-party members",
  "The state legislature bans all out-party members from competing in state and local campaigns in response to recent controversies." = "Bans all\nout-party members",
  "Your local governing board is evenly split between Democrats and Republicans." = "Local board\nevenly divided",
  "In-party takes control of your local governing board." = "In-party takes\nlocal board",
  "Out-party takes control of your local governing board." = "Out-party takes\nlocal board",
  "In-party wins the Presidency and takes control of Congress." = "In-party wins Presidency\n& Congress",
  "In-party wins the Presidency but Congress is divided." = "In-party wins Presidency,\ndivided Congress",
  "Out-party wins the Presidency but Congress is divided." = "Out-party wins Presidency,\ndivided Congress",
  "Out-party wins the Presidency and takes control of Congress." = "Out-party wins\nPresidency & Congress",
  "The state legislature leaves your pay unchanged." = "Pay unchanged",
  "The state legislature cuts your pay by 20%." = "Pay cut by 20%",
  "The state legislature cuts your pay by 10%." = "Pay cut by 10%",
  "The state legislature cuts your pay by 5%." = "Pay cut by 5%",
  "The state legislature raises your pay by 5%." = "Pay raised by 5%",
  "The state legislature raises your pay by 10%." = "Pay raised by 10%",
  "The state legislature raises your pay by 20%." = "Pay raised by 20%"
))

mm_resultsPID_df$Level <- factor(mm_resultsPID_df$Level, levels = names(new_labels), labels = new_labels)
table(mm_resultsPID_df$Level, useNA = "always")


ggplot(mm_resultsPID_df, aes(x = Level, y = Estimate, color = Party)) +
  geom_point(position = position_dodge(width = 0.75), size = 1.5) +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), 
                position = position_dodge(width = 0.75), 
                width = 0.1,
                linewidth = 0.5) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  coord_flip() +
  facet_wrap(~Attribute, scales = "free_y", ncol = 1) +
  labs(x = "", y = "Marginal Mean", title = "") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 12, color = "black"),          
    axis.text.y = element_text(size = 12, color = "black"),          
    axis.title.x = element_blank(),
    axis.title.y = element_text(size = 14, color = "black"),         
    strip.text = element_text(size = 13, angle = 0, color = "black"), 
    panel.grid.major = element_blank(),             
    panel.grid.minor = element_blank(),             
    panel.background = element_rect(fill = "white", colour = "white"),
    axis.line.y = element_line(color = "black"),
    axis.ticks.y = element_line(color = "black"),
    axis.line = element_line(colour = "black"),
    legend.text = element_text(size = 12, color = "black"),          
    legend.title = element_text(size = 12, color = "black"),         
    legend.position = "bottom"                      # Move legend to bottom
  ) +
  scale_color_manual(values = c("Democrat" = "blue", "Republican" = "red")) +
  scale_x_discrete(labels = new_labels)  

ggsave("Figure2.pdf", width = 6, height = 11, units = "in", dpi = 500)


## Education 

long$CollegeEducation <- as.factor(long$CollegeDegree)

mm_resultsCollege <- cj(long, chosen ~ Election_Legislation + `2024_Local_Election_Outcome` + `2024_Presidential_Election_Outcome` + Salary_Change, estimate = "mm", id = ~ ID, by = ~ CollegeEducation)
